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ABSTRACT 

This paper presents a detailed statistical determination of the equatorial rotation rates of 
classical Be stars. The rapid rotation of Be stars is likely to be linked to the ejection of gas that 
forms dense circumstcllar disks. The physical origins of these disks are not understood, though it 
is generally believed that the ability to spin up matter into a Keplerian disk depends on how close 
the stellar rotation speed is to the critical speed at which the centrifugal force cancels gravity. 
There has been recent disagreement between the traditional idea that Be stars rotate between 
50% and 80% of their critical speeds and new ideas (inspired by the tendency for gravity darkening 
to mask rapid rotation at the equator) that their rotation may be very nearly critical. This paper 
utilizes Monte Carlo forward modeling to simulate distributions of the projected rotation speed 
(usini), taking into account gravity darkening, limb darkening, and observational uncertainties. 
A chi-squared minimization procedure was used to find the distribution parameters that best 
reproduce observed vsmi distributions from R. Yudin's database. Early-type (07e-B2e) Be stars 
were found to exhibit a roughly uniform spread of intrinsic rotation speed that extends from 40- 
60% up to 100% of critical. Late-type (B3e-A0e) Be stars exhibit progressively narrower ranges 
of rotation speed as the effective temperature decreases; the lower limit rises to reach critical 
rotation for the coolest Be stars. The derived lower limits on equatorial rotation speed represent 
conservative threshold rotation rates for the onset of the Be phenomenon. The significantly 
subcritical speeds found for early-type Be stars represent strong constraints on physical models 
of angular momentum deposition in Be star disks. 

Subject headings: circumstellar matter — stars: atmospheres — stars: early-type — stars: 
emission-line, Be — stars: fundamental parameters — stars: rotation 



1. Introduction 

Be stars are rapidly rotating, non-supergiant B-type stars that exhibit, or have exhibited in the past, 
emission in their hydrogen Balmer lines. The observed properties of Be stars are consistent with the coexis- 
tence of a dense circumstellar disk (flattened in the plane perpendicular to the rotation axis) and a variable 
stellar wind (Struve 1931; Doazan 1982; Slettebak 1988; Prinja 1989; Porter & Rivinius 2003). The gas in 
the so-called "decretion disk" is traditionally believed to be ejected from the star and not accreted from an 
external source (see, however, Harmanec et al. 2002; Abt 2004). Although there is increasing evidence that 
the disk gas is in Keplerian orbit (e.g., Hanuschik 1996; Hummel & Vrancken 2000), there is a great deal 
of evidence that Be-star photospheres are rotating too slowly to propel any atmospheric material into orbit. 
Typical observationally determined values of the ratio of equatorial rotation speed to the critical rotation 
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speed V cr it (at which gravity is balanced by outward centrifugal forces) range between 0.5 and 0.8 (Slettebak 
1982; Porter 1996; Yudin 2001). If this is the case, then any theoretical model for the origin of Be-star disks 
would require a substantial increase in angular momentum between the photosphere and the inner edge of 
the disk. 

Recently, the idea that Be stars are rotating with significantly subcritical rotation speeds has been 
called into question. The primary observational diagnostic of hot-star rotation is the Doppler broadening of 
photosphcric absorption lines, first elucidated by Abney (1877). Rotational broadening provides a surface- 
weighted measure of the product of the equatorial rotation speed V cq and sinz, where i is the inclination 
angle to the observer. Traditional means of determining V4 q sinz from line profiles (e.g., Tassoul 1978; Gray 
1992) often assume that the star is spherical. However, rapidly rotating O and B stars tend to become 
centrifugally distorted into oblate shapes and thus undergo "gravity darkening" (i.e., a redistribution of 
radiative flux in proportion to the centrifugally modified gravity; von Zeipel 1924). The equators of such 
distorted stars become dimmer and cooler than their poles, and thus the most rapidly rotating regions of the 
stellar surface are weighted less strongly in the resulting star-averaged absorption profiles. This tendency 
for gravity-darkened stars to exhibit narrower profiles than would be the case for spherical stars — and thus 
lower computed values of V4 q sini — has been known for more than a half century (e.g., Slettebak 1949; 
Stoeckley 1968; Hardorp & Strittmatter 1968; Walker et al. 1979; Collins & Truax 1995) and has been 
recently highlighted as a potential bias in statistical samples of Be star rotation rates (Zorec et al. 2003; 
Townsend et al. 2004; Cohen et al. 2005; Frcmat et al. 2005). 

A physical understanding of the Be phenomenon hinges on how close the stars are rotating to their 
critical speeds. If V cq is within one or two sound speeds of V a i t (which would imply V eq /V a i t > 0.95), there 
are many possible weak processes that could easily propel gas into orbit (see, e.g., Owocki 2005). When the 
above ratio falls below ~0.9, though, the increased amount of energy and angular momentum addition that 
would be needed to spin up material into a Keplerian disk is large enough to greatly restrict the number 
and type of potential sources. Townsend et al. (2004) suggested that gravity darkening effects could be 
strong enough to make a distribution of nearly critical rotation speeds appear to be shifted down to values 
of 0.5-0.8Krit if the line profiles were interpreted as if the stars were spherical. The inclusion of gravity 
darkening, however, tends to complicate the analysis to the extent that a unique determination of V eq sin i 
from a single measured line width (for any individual star) does not seem to be possible. 

This paper attempts to disentangle the above effects by using Monte Carlo forward modeling to produce 
a large number of trial probability distributions of V cq . Each distribution is processed, assuming random 
inclination angles (and with inclination-dependent line narrowing due to gravity darkening), to simulate an 
observed statistical sample of line widths. The most likely intrinsic distribution of Be-star rotation speeds 
is thus determined by searching for the models with the minimum x 2 differences between the simulated and 
observed line width distributions. The derived distributions of V eq , as a function of spectral type, yield 
important empirical constraints on the threshold rotation speeds for the occurrence of the Be phenomenon. 
This forward-modeling method is less ambiguous than the more common inverse technique of using simple 
geometric transformations to convert an observed distribution of V eq sin i values into either a distribution of 
intrinsic rotation speeds or a mean value for V cq . 

Although this kind of analysis has a long history (e.g., Chandrasekhar & Munch 1950; Stoeckley 1968; 
Lucy 1974; Balona 1975; Porter 1996; Clark & Steele 2000; Chauville et al. 2001), the present work contains 
several novel features that help to increase the overall level of confidence in the results. First, the number 
of observed stars — from the published database of Yudin (2001) — is now large enough to be able to use 
the detailed shapes of the number distributions as constraints rather than just their low-order moments. 
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Second, the effects of gravity darkening are included in the most "conservative" manner possible, thus 
taking into account the heterogeneous origins of the V eq sini entries in the database (i.e., gravity darkening 
was considered in the calculation of some V eq sini values, but not others). The resulting subcritical values of 
V cq are thus designed to be safe upper limits, and the actual rotation speeds may be even lower if the modeled 
gravity darkening effects were overestimated. Third, the derived rotation speeds are used as inputs to an 
independent statistical simulation of visible polarization measurements of Be-star disks. The good agreement 
between the shapes of the observed and simulated distributions of polarization is a useful validation of the 
derived range of V eq values. 

The remainder of this paper is organized as follows. § 2 presents a summary of the Yudin (2001) Be-star 
database and a description of the adopted fundamental stellar parameters that were used to compute V C Ht 
and other physical quantities. § 3 describes the Monte Carlo forward modeling procedure that was used to 
simulate statistical distributions of Be stars, and also gives the resulting best-fit ranges of equatorial rotation 
speed. The possibility that nonstandard gravity darkening exponents may apply to Be stars is investigated in 
§ 4, and a simulation of linear polarization values for the derived distribution of rotation speeds is presented 
in § 5. Additional pieces of evidence in favor of the results derived in § 3 (essentially that V cq ^ V CI { t for 
early-type Be stars) are laid out in § 6. A summary of the major results of this paper, together with a 
discussion of the implications for theories of the Be phenomenon, is given in § 7. 



The Yudin (2001) database of early-type emission-line stars contains 627 objects with MK spectral types 
between 07. 5e and B9/A0e and luminosity classes between II/III and V. Care was taken to exclude Of, B[e], 
and Herbig Ae/Be stars from this database, thus making it the largest sample of "classical" Be stars yet 
assembled. There are 462 stars in the catalog 1 with nonzero values for the projected rotation speed wsini, 
and this subsample contains the primary observational data to be compared with the Monte Carlo model 
predictions in § 3. (The lower-case v is used here only in combination with sini to denote the convolved 
quantity derived empirically from line widths. The nomenclature V cq sini is used below only for the product 
of two known quantities.) 

In order to remove potential biases arising from the substantial variation of stellar parameters from the 
late-0 to early- A spectral ranges, the observed values of vsini for each star should be normalized by the 
star's critical rotation speed. Fundamental parameters (e.g., mass M* and polar radius R p ) for each star are 
needed to compute the critical rotation speed, which is defined for a rigidly rotating Roche-model star as 



with G being the Newtonian gravitation constant (see, e.g., Jeans 1928; Collins 1963; Tassoul 1978). The 
above expression is consistent with the existence of continuum radiation pressure as long as von Zcipcl (1924) 
gravity darkening applies and the continuum Eddington factor T is less than ^0.5 (Glatzel 1998; Maeder & 



1 Note that Yudin (2001) states that there are 463 stars that have nonzero values of the projected rotation speed, but the 
online version of the database (VizieR catalog J/A+A/368/912) appears to contain only 462 nonzero values of vsini. This 
discrepancy is unimportant for any of the statistical results of Yudin (2001) or this paper. 



2. 



The Observational Database 




(1) 
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Meynet 2000a). For completeness, the Eddington factor is given by 

1 " 4,cGM, ~ 2 - ? 1U (M*/M ) ' (2j 

where c is the speed of light in vacuum, a e is the Thomson scattering opacity, and M Q and L Q are the Sun's 
mass and bolomctric luminosity. The numerical factor above was computed from a standard solar abundance 
mixture (X = 0.73, Y = 0.24). For B-type stars, T is typically much smaller than 1, thus justifying the 
choice of solution branch to the radial force balance equation (at the surface of the critically rotating star) 
that is implied by equation (1). The largest value of T in Yudin's (2001) database is 0.23 (for the 07.5 III 
star 68 Cyg), and there are only 3 other stars out of 627 that have T > 0.1. 

The spectral types and luminosity classes listed in Yudin's (2001) catalog were used to compute bolo- 
mctric luminosities L* and mean effective temperatures T e g by using the statistical relations of de Jager & 
Nieuwenhuijzen (1987). Their mean tabulated values were derived from a set of 199 high-precision determi- 
nations of L* and 268 high-precision determinations of T e g across the Hertzsprung-Russell (H-R) diagram. 
For uncertain spectral types and luminosity classes that were listed by Yudin (2001) using two possible values 
(e.g., "B9/A0" or "III/IV"), the luminosities and temperatures were computed for each value then averaged 
logarithmically. Stars without a listed luminosity class were assumed to be main sequence (class V) objects. 
Although more recent calibrations of L* and T e g exist for the earliest-type stars (see, e.g., Garcia & Bianchi 
2004; Crowther 2005), the de Jager & Nieuwenhuijzen (1987) tables remain the most complete and cohesive 
set of correspondences that is separated by luminosity class. 

Once L* and T e g were determined for each star, the stellar mass M* was computed by interpolating 
between the evolutionary tracks published by Claret (2004). First, the abscissa in the H-R Diagram was 
transformed from T e g to the scaled variable S = T c g/L® A6 , which skews the main sequence to be roughly 
vertical. Each of the 30 evolutionary tracks (spanning initial masses between 0.8 and 126 M ) was searched 
for the point where 5 matched that of the star in question. The luminosities and masses at these points were 
saved into one-dimensional arrays, thus effectively giving a tabulation of M* as a function of L t . The actual 
stellar mass was then computed by linear interpolation using the empirically determined value of L*. The 
resulting mass-luminosity relationships for luminosity classes III, IV, and V were fit with simple quadratic 
functions (in logarithm space) and are given here for comparison to other calibrations: 

y = 0.2084 + 0.0634a; + 0.0325a; 2 (III) (3) 

y = 0.2055 + 0.0746a; + 0.0322a; 2 (IV) (4) 

y = 0.1997 + 0.0844a; + 0.0312a; 2 (V) (5) 

where y = log 10 (M»/M©) and x = log 10 (i*/L Q ). These fits apply only to the late-0 to late-B range 
represented in Yudin's (2001) database (i.e., masses between 2.5 and 38 M Q ) and the fits are accurate to 
within 5% in the mass over this range. 

Figure 1 shows the spectral type dependence of V cr ; t as computed from equation (1) for luminosity 
classes III, IV, and V. Corresponding values of the critical rotation speed given by Porter (1996) and Yudin 
(2001) are also shown. For main sequence stars, the values computed for this paper are in good agreement 
both with these other plotted values and with the often-cited tables of Schmidt-Kaler (1982), Underhill 
(1982), Harmanec (1988), and Andersen (1991). For luminosity classes III and IV, the computed V cr i t values 
are systematically larger than those of Yudin (2001). This is the result of the trend that the luminosities 
of de Jager & Nieuwenhuijzen (1987) tend to be on the low side when compared with other calibrations 
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for giants and subgiants. For constant T c g, the stellar radius computed with a lower luminosity would be 
smaller, and thus V cr it would be larger. 

Figure 1 also displays the V CI it values computed by Chauville et al. (2001) for each of the 116 Be stars 
in their published database. These values are plotted with the same symbols for all luminosity classes, and 
it can be seen clearly that they are about 20% lower than the values derived in this work and on average 
they are even slightly lower than Yudin's (2001) values. The reasons for this systematic discrepancy are not 
clear, although it seems to be at the root of the significant difference between the derived values of the mean 
ratio usini/V^ r it of Yudin (2001) and Chauville et al. (2001); see below. 

There are two main pieces of evidence that support the comparatively large values for V cr it derived in 
this paper: 

1. Recent model atmosphere based determinations of B-star fundamental parameters (Fitzpatrick & 
Massa 2005) give critical rotation speeds that agree well with the curves shown in Figure 1. For 
the 5 luminosity class III stars (ranging from B3 to AO) in the Fitzpatrick & Massa (2005) database, 
the computed values of V C rit were all slightly larger than the values given in Figure 1. For the 20 (17) 
stars of class IV (V), the computed y cr i t values fell roughly equally above and below the respective 
curves in Figure 1. 

2. When modeling the properties of Be stars above the main sequence (i.e., luminosity classes III and 
IV) , it is probably safest to choose the set of critical rotation speeds that cleaves closest to the main 
sequence values. It has been known for some time that rapid rotation and gravity darkening can make 
a main sequence star appear to be up to a full magnitude brighter than than its surface-averaged 
luminosity would indicate (e.g., Maeder & Peytremann 1970; Collins & Sonneborn 1977; Collins et al. 
1991). Despite some evidence that the Be phenomenon occurs only during relatively late evolutionary 
phases, it is possible that the true distribution of Be stars on the H-R diagram remains closer to the 
main sequence than is implied by the inferred fractions of luminosity class III and IV stars (see also 
Fabregat & Torrejon 2000). Therefore, it may err on the side of caution to use the present calibration, 
which yields relatively high values of V cr i t for Be stars classified as III and IV. 

More work needs to be done to more accurately pin down the fundamental parameters of Be stars. By 
default, the statistical models presented below mainly use the high-end y cr it values derived for this paper 
(sometimes denoted as "Cranmer"), but a parallel analysis is also done using the low-end Chauville et al. 
(2001) values. 

The 462 stars with measured projected rotation speeds were normalized by dividing Yudin's (2001) 
tabulated vs'mi by the computed values of V cr it for each star. The mean value of the ratio v smi/V CI it was 
found to be 0.485, which is slightly smaller than the value of 0.50 found by Yudin (2001) using lower critical 
rotation speeds. Also, the latter value was computed as the ratio of the mean vsini to the mean V a { t , 
not as the mean value of the set of ratios. The dimensionless standard deviation, skewness, and kurtosis of 
the distribution of ratios (as defined by Press et al. 1992) were found to be 0.167, —0.00106, and —0.432, 
respectively. The significantly negative kurtosis denotes a flat-topped plateau-like distribution that seems to 
occur because the full sample of stars is made up of subpopulations with a range of mean values that sweep 
across the central peak of the full distribution (see below). This result highlights the usefulness of using the 
ratio vsinz/T4 r it rather than just vsini itself, because Yudin (2001) found that the observed distribution of 
v sin i values was well- represented by a normal distribution. Thus the "normality" of the v sin i distribution 
could be related more to the large spread in V CI i t across the B-type spectral range than to the true distribution 
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of relative rotation rates. 

To facilitate the study of the onset of the Be phenomenon as a function of spectral type, the 462 stars 
in the database were divided into five subpopulations, two containing 93 stars and three containing 92 stars. 
Table 1 contains a summary of the properties of these subsamplcs, which were generated by binning stars 
of similar T eT \. From hot to cool effective temperatures, the groups are labeled very early (VE), early (E), 
medium (M), late (L), and very late (VL). The exact values of T e g used for the dividing lines between the 
groups were chosen in order to give all five subpopulations a roughly equal number of stars. Table 1 also lists 
the mean, standard deviation, skewness, and kurtosis of the vsmi/V cr it distributions for each subpopulation. 
For further statistical analysis as a function of spectral type and luminosity class, see Yudin (2001). 

Figure 2 displays the number distribution of usini/V^ r i t for all 462 Be stars, as well as the distributions 
for the five subpopulations defined in Table 1, each normalized to an integral of unity (i.e., expressed as 
binned probability densities). As noted by Yudin (2001), the mean values of v smi/V cr it have a systematic 
dependence on T c g, increasing monotonically from ~0.4 for the hottest (VE) stars to ~0.6 for the coolest 
(VL) stars. There is no clear spectral type dependence on the standard deviations (or higher moments) of the 
subpopulation distributions. The distributions are plotted by dividing the region between vsmi/V CI i t = 
and 1 into 11 equally spaced bins; this provides a good balance between resolution and statistics. Note that 
the smallest and largest computed values of vsini/V cr i t for the full database are 0.035 and 0.912. (The fact 
that there are no stars with measured v sin i that exceed V CI i t tends to support the conclusion reached below 
that Be stars are not all rotating critically; see also § 4). 

As an alternate derivation of the normalized projected rotation speeds, the 462 wsini values were also 
divided by the lower V cr i t values indicated by Chauville et al. (2001). A least-squares power-law fit was found 
for the spectral type dependence of the critical rotation speeds from that paper: V CI it = (450 km s _1 )/s ' 775 , 
where s is a continuous variable that straightforwardly denotes the spectral subtype (for 05, B0, and B9, 
s = 0.5, 1.0, and 1.9, respectively). This fit does not reflect the spread in values from the range of luminosity 
classes (see Figure 1) but it accurately models the trend for the Chauville et al. critical speeds to be lower 
than those derived by others. Table 1 gives the mean, standard deviation, skewness, and kurtosis of the 
t)sini/T4 r it distributions for each subpopulation using this alternate calibration for V crl f For the earliest 
spectral bin, VE, the mean ratio of 0.49 is only slightly higher than the value of 0.41 that corresponds to 
the V cr it calibration derived above. The differences between the two calibrations grow steadily larger for 
the cooler spectral bins, leading to mean ratios of either 0.90 (Chauville et al.) or 0.59 (Cranmer) for the 
VL subpopulation, depending on V^rit- For the entire database, the mean ratio using the Chauville et al. 
calibration for V CI n is 0.627; this is very close to the value of 0.65 given by Chauville et al. (2001) for their 
sample of 116 Be stars. This value is almost 30% higher than the mean ratio of 0.485 that was computed 
above, and the difference is due solely to the differences in stellar mass and radius that go into the calculation 
of the critical rotation speed. 



3. Forward Modeling of Distributions of the Projected Rotation Speed 

In this section the procedure for simulating theoretical distributions of v ami/V cr it for a random sample 
of Be stars is described. The simulation of wsini measurements for individual gravity darkened stars is 
outlined in § 3.1 (see also Collins 1974; Collins & Truax 1995; Townsend et al. 2004). The construction 
of theoretical probability distributions for a large number of stars, together with the comparison with the 
observed distributions, is described in § 3.2. 
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3.1. Modeling Individual Stars 

For a given set of fundamental stellar parameters (mass, luminosity, and polar radius), an equatorial 
rotation rate (V eq ), and an observational inclination angle (i), it is possible to compute unique absorption 
line profile shapes that can be used as direct measures of the quantity v sini/V cr it- Figure 3 illustrates the 
effect of oblatcness and von Zeipel (1924) gravity darkening on the absorption line shape for a representative 
late-B main sequence star (see also Cohen et al. 2005). As expected, gravity darkening preferentially weights 
the more slowly rotating polar regions and thus produces a narrower line profile than that computed without 
gravity darkening. To produce the images and profiles in Figure 3, the computer code from Cranmer & 
Owocki (1995) was used to generate high-resolution models of both a spherical star and an oblate Roche- 
model star with the same polar radius and equatorial rotation velocity. The stars are observed from a large 
distance (1000 R p ) at a given inclination angle, and the same limb darkening law is applied in both cases 
(see below). Each pixel of the image has a projected linear dimension of 0.008 R p in the plane perpendicular 
to the line of sight that crosses through the center of the star. 

The rotationally broadened profiles, here assumed to be the commonly used He I A4471 absorption lines, 
were simulated for the two cases shown in Figure 3 by the following procedure. The residual flux is defined 
as the ratio of the flux in the spectral line to that of the surrounding continuum, 

n = J dxdyl„, h {x,y) 
J dxdy I u . c (x,y) 

where x and y are spatial coordinates that specify a grid of rays parallel to the observer's line of sight. For 
the subset of these rays that intercept the stellar surface, Cranmer & Owocki (1995) and Cranmer (1996) 
described how to convert these coordinates into star-centered spherical coordinates r and 9 (ignoring the 
azimuthal longitude (f> because the stars are assumed to be axially symmetric). The colatitude determines 
the local magnitude of the centrifugally modified effective gravity, and also — for the oblate gravity-darkened 
case — the local effective temperature. 

For simplicity, the specific intensity of the continuum is modeled as the product of a Planck function at 
the local effective temperature and a linear limb darkening function, i.e., 

5 

h,c = -i-B u (T cB ) (1 - u c + u c fi) , (7) 

47T 

where u c is the continuum limb darkening constant and [i is the cosine of the angle between the line of 
sight and the normal to the stellar surface. Below, we adopt u c = 0.4 from both a classical treatment of 
limb darkening for a gray temperature distribution at a continuum wavelength of 4471 A (see eq. 10-20 of 
Mihalas 1978) and a numerical study of limb darkening for rapidly rotating B stars (Table 4 of Collins & 
Truax 1995). The specific intensity in the spectral line /„! is modeled as the product of I Vj c and a residual 
intensity /„ that is computed from Milnc-Eddington theory for a pure absorption line, with 

l-U c + U c ^/{l + T] v ) 

Jv — : j (°) 

1 - U c + U c [l 

(Mihalas 1978; Collins 1989). The dimensionless line absorption coefficient t] v is modeled as a Voigt profile 
that is thermally broadened and Dopplcr shifted by the stellar rotation, 

Vv = voH(a,0 , (9) 

and 

^-^(I + Vlos/c) 
5 V th v /c [W> 
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w = 2u c is v th r°° 

c[(l - u c )V3 + u c ] 7o 



(11) 



where £f(a,£) is the standard Voigt function, i^o is the rest-frame frequency of the line transition, Vlos is 
the projected component of the rotation velocity along the line of sight for this ray, and Vth is the thermal 
speed of helium atoms corresponding to the local effective temperature. The dimensionless strength of the 
Voigt wings is assumed here to be given by a constant value of a — 0.08 based on standard Stark broadening 
for the He I A4471 line at a mean B-star temperature of 2 x 10 4 K (e.g., Griem et al. 1962; Barnard et al. 
1969; Leckrone 1971). 

The strength of the absorption line for each ray that intercepts the stellar surface is parameterized by 
the dimensionless line-center absorption coefficient t/q. This quantity is determined empirically from a grid 
of detailed spectral synthesis calculations of the He I A4471 equivalent width from a collection of early- 
type atmosphere models (Gonzalez Delgado & Leitherer 1999). The effective gravity and temperature at 
the stellar surface for each ray is used to interpolate the local equivalent width from Table 6 of Gonzalez 
Delgado & Leitherer (1999). The conversion between r]o and equivalent width W is given by the theoretical 
curve of growth, 

VoH(a,Q 
_l + rioH(a,Z)_ 

where W above is in frequency units (Mihalas 1978). This relation has been computed numerically and 
tabulated for the adopted values of a and u c . For B stars, the equivalent widths for individual rays span the 
range between 0.1 and 1.5 A, corresponding roughly to rjo values between 30 and 3000. 

The procedure outlined in equations (6)— (11) is certainly more simplistic and approximate than per- 
forming full model atmosphere calculations of the relevant line profiles (as was done by, e.g., Townsend et al. 
2004). However, the relative ease of computing reasonably accurate profile shapes by the above method allows 
very fine grids of wavelengths, inclination angles, and rotation speeds to be computed without prohibitive 
computational expense. These grids make possible the detailed statistical studies described below. 

Spectral lines were computed for representative main sequence stars of spectral types B0, B2, B3, B5, 
and B9; i.e., for stars roughly at the centers of the five subpopulation bins defined in § 2. For each spectral 
type, line profiles for the spherical and the oblate gravity darkened cases were computed on a fine grid of 
200 wavelength points, and for 100 inclinations (i between 0° and 90°) and 100 equatorial rotation speeds 
(Veq/V cr it between and 1). A single measure of the line width was determined from the full width at 
half maximum (FWHM) of the numerically computed profiles. The key parameter to be used below is the 
dimensionless ratio D, which is defined as the FWHM line width for a model star computed with oblateness 
and gravity darkening divided by the corresponding FWHM for the spherical model. Thus, D provides 
an indication of the relative change in the absorption line width that is produced by the effect of gravity 
darkening. Figure 4 shows D as a function of both Ksq/Krit and i for the earliest and latest (B0 and B9) 
spectral types. The information contained in Figure 4 is essentially equivalent to that shown in, e.g., Figure 
4 of Stoeckley (1968), Figure 6 of Collins & Truax (1995), Figure 1 of Townsend et al. (2004), and Figure 7 of 
Fremat et al. (2005), but for different spectral types. These other figures have tended to plot the computed 
line widths as functions of projected rotation speeds, which highlights the ambiguity involved with attempting 
to determine the product V4 q sini from from an "observed" wsini. The forward modeling procedure outlined 
in this paper is essentially free of such ambiguity. 

The overall impact of the line narrowing due to gravity darkening can be understood better by converting 
from D to the similarly defined "velocity deficiency" quantity of Townsend et al. (2004), 

6V = (1 - D)V eq smi/V CTit , (12) 

specifically for their fiducial case of edge-on inclination (i = 90°) and T4q/Krit = 0.95. For the five spectral 



- 9 - 



subtypes, there is a gradual increase of the effect as one goes from BO (SV = 12.5%) to B9 (SV — 26.0%), 
with intermediate values of 16.0%, 17.5%, and 21.0% for B2, B3, and B5, respectively. These values compare 
favorably with those given in Table 2 of Townsend et al. (2004), which helps to validate the use of the simpler 
line synthesis technique described above. 

There are noticeable differences between the shapes of the contours in the B0 and B9 cases shown in 
Figure 4. The intermediate (B2, B3, B5) spectral types have contours in D that change in shape gradually 
between the two plotted extreme cases. These differences arise because the He I A4471 equivalent widths, as 
interpolated from the modeled spectra of Gonzalez Delgado & Leitherer (1999), depend on latitude differently 
for different ranges of effective gravity and temperature. An early (B0) model rotating at 99% of critical 
exhibits a maximum in W at mid-latitudes of about 0.9 A, and minima at the poles and equator of ~0.75 
and ^0.3 A, respectively. A similar late-type (B9) model exhibits a simpler monotonic decrease in W from 
pole (^0.35 A) to equator (~0.1 A). These differences result in different spectral line shapes and different 
intensity weightings over the oblate surfaces. 

The procedure to simulate a measurement of vsmi/V cr i t for an individual star is summarized by the 
following relation: 



The choice of a range of values for V eq /V cr i t is discussed below in § 3.2, and the inclination angle i is defined 
formally as cos _1 /ii, where /i^ is chosen randomly from a uniform probability distribution between and 
1 (see also Chandrasekhar & Munch 1950). The gravity darkening factor D is interpolated from the grids 
of values that were used to generate Figure 4. The factor ( above is a simulated observational uncertainty, 
which is sampled from a normal random distribution having a mean of zero and a standard deviation of 
(i.e., 68% of the time, £ falls between — <r<; and +<7<j). The "mean uncertainty level" <7<j is kept constant for 
any given sample of model stars, though the effects of varying this parameter between and 0.3 are explored 
below (see also Balona 1975). Yudin (2001) discussed the determination of standardized error bounds for the 
462 observed v smi values in the database, and found typical relative uncertainties of order 10%. Uncertainty 
levels up to 30% in the ratio vsmi/V CI i t may be reasonable to assume, since neither the numerator nor the 
denominator are known precisely. 

In order to most stringently test the hypothesis that Be stars are rotating nearly critically, it would be 
desirable to make the assumptions that are most favorable to that hypothesis (i.e., those that tend to give 
largest derived values of Kq/Krit)- Thus, if the derived rotation speeds are still substantially below critical 
even when those assumptions are made, the hypothesis of critical or nearly critical rotation can be ruled out 
at a high level of confidence. Two such assumptions regarding the use of D in equation (13) are adopted 
here: 

1. Despite the variations of D as a function of spectral type, in the statistical models below we apply 
the grid of D values computed for the B9 case to all of the spectral-type subpopulations. This is 
done because the B9 case shows the strongest line-narrowing effect due to gravity darkening (i.e., 
the lowest values of D for rapid rotation and edge-on inclination). Applying this lower- limit case for 
D in equation (13) leads to a possible overestimate of the narrowing trend and thus a systematic 
underestimate of vsmi/Vc^ for a specific assumed V^ q . When compared with observed values of 
vsmi/Vcrit, then, the resulting derived V cq would end up being an overestimate. 




(13) 



2. The fact that D is included at all in the equation above implies that the reported measurements 
of vsmi/V CT n have not taken gravity darkening into account. However, several of the sources of 
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observational v sin i data used by Yudin (2001) certainly did their best to account for gravity darkening, 
and thus they reported processed estimates of V eq smi rather than raw measurements of vsini. If this 
is the case for a substantial fraction of the observations, then the use of equation (13) would lead to a 
systematic underestimation of wsini/y cr i t . This effect works in the same sense as the overestimation 
of the strength of line narrowing from using the B9 models for all spectral types. 

Thus, the present models stand at one end of a continuum of modeling choices (i.e., "strong" gravity darkening 
effects), and the assumption of no gravity darkening (i.e., D = 1) would stand at the opposite end. The true 
rotation speeds of the Be stars should fall somewhere between the values derived for these limiting cases. 



There have been numerous attempts to deconvolve intrinsic statistical distributions of stellar rotation 
speeds from the observed distributions of vsini values. Analytic studies include Chandrasekhar & Munch 
(1950), Bernacca (1970, 1972), Lucy (1974), and Balona (1975). More recent efforts to determine the 
statistical distributions numerically — rather than just estimate mean values or assume that all stars have 
the same rotation rate — include Wolff et al. (1982), Chen & Huang (1987), Porter (1996), and Brown & 
Verschueren (1997). Generally, these models tend to be "inverse" determinations that begin with an observed 
distribution function of projected rotation speeds $ (vsini) and work iteratively towards a consistent form 
for the true distribution function of intrinsic rotation speeds ty(V eq ). This technique becomes potentially 
ambiguous, though, when the line-narrowing effects of gravity darkening are taken into account; i.e., how 
are we to be sure that the derived solution for ty(V eq ) is unique when an observed vsini cannot be mapped 
identically to a single value of the product I4 q sini? 

Although the present forward-modeling method docs not produce completely unique solutions, it treats 
all possible distributions on the same footing and thus is not in danger of missing potential solutions. After 
choosing a parameterized functional form for the distribution of rotation speeds ^(V cq ), the parameters are 
varied by producing a multidimensional grid of trial distributions having the full range of combinations of 
parameters. (All velocities here are assumed to be in units of V^ r it-) Each trial distribution is converted into 
an observed distribution $(vsini) by using equation (13), and each of these is compared to the observed 
distributions for the five spectral-type subsamples shown in Figure 2. The best fits are judged with the 
X 2 diagnostic appropriate for comparing a coarsely sampled observed distribution with a known model 
distribution (see below). 

Earlier studies have used various functional forms for ty(y eq ) such as Dirac delta functions, Gaussians, 
and uniform (i.e., flat-topped) distributions. After some experimentation, it was found that a truncated 
linear function of the form 



best balanced the demands of versatility and simplicity. (From their appearance when plotted, these dis- 
tributions can also be called "trapezoidal.") The three free parameters of this family of functions are the 
minimum and maximum truncation values (Vmin and V max ) and the slope to. For any choice of these pa- 
rameters, the constant b is determined automatically by assuring that ^(V eq ) is normalized to unity upon 
integration over all V eq . The minimum and maximum possible slopes are defined by the distribution being 
"right triangular;" i.e., the minimum slope occurs when ^(Vmin) is the global maximum and ^(Vmax) = 0. 



3.2. 



Monte Carlo Distributions and Results 




Knin < V cq < V Tj 

otherwise, 



max 



(14) 



- 11 - 

The maximum slope occurs when ^(Knax) is the global maximum and ^(T^min) = 0. Thus it is convenient 
to define a dimensionless parameter S that ranges between —1 and +1, with the slope ranging from its 
minimum to its maximum value over this range, and 

_ 25 



Kiax Vmin 

The parameter S tends to give a better qualitative impression of the shape of the function than does m. 

For each choice of the above parameters, the resulting distribution function ^ then needs to be converted 
into a distribution function $ for the observed values of v sini/V^rit- Rather than simulating large numbers 
of stars for each point in the three-dimensional grid of parameter space (V m i n , V max , S), a series of Monte 
Carlo simulations was performed for a one-dimensional grid of delta-function distributions, 

W=^-^). (17) 

with 80 equally spaced values of Vj/V cr i t ranging between and 1. These distributions then served as basis 
functions used to build up any desired arbitrary statistical distribution. Each Monte Carlo simulation used 
10 5 model stars with random inclinations i and observational uncertainties ( as described in § 3.1, and 
identical values of V cq — Vj. Equation (13) was used to simulate individual values of vsmi/V CI it, and these 
values were summed into the same 11 bins as were used in Figure 2 for comparison with observations. This 
process yielded a a one-dimensional grid of <Fj(usini) functions corresponding to the delta- function tyj grid. 
Finally, then, the <& function corresponding to a given trapezoidal distribution (eq. [14]) was computed by 
taking the linear combination 

$(wsini) oc ^ ®j(vsmi) *(Vj) , (18) 

3 

and renormalizing to an integral of unity. In other words, the trapezoidal distribution is used as a weighting 
function to determine what fraction of each of the 80 $j functions are summed together to form the complete 
distribution $. 

Before presenting the results of the above process, it is instructive to show the $j(usini) distributions 
that correspond to Vj = V CT it. These represent the predicted distributions of projected rotation speed 
assuming that all Be stars rotate critically. Figure 5 shows a series of these distributions that were computed 
for different choices for the mean uncertainty level o~q. For o~q = 0, the maximum modeled value of v sini/V^nt 
is ^0.68, corresponding essentially to the maximum value of Dsini along the right edge of Figure 4b. (If 
gravity darkening had not been taken into account, the maximum modeled value of vsini/V CI i t for the 
<7£ = model would have been 1.) Note the secondary peak around vsmi/V cr i t w 0.53; this corresponds 
to the small local maximum in D for i w 40° and V eq /V CI it = 1. Predictably, when <tq is increased, the 
probability distributions become smeared out into increasingly Gaussian shapes. 

Figure 5 also shows the observed $ distribution of the earliest-type (VE) subpopulation of Be stars, 
normalized by the low choice of V CT it from the fit to the Chauville et al. (2001) values. (The distribution is 
shifted to the right compared to the one shown in Figure 2b.) Thus, the observed distribution is plotted 
with assumptions that push the peak to high values of v sini/V^ r it and the modeled distributions are plotted 
with the "strong gravity darkening" assumptions that push the peaks to low values. From the fact that 
there is still a 15-20% mismatch between the observed and modeled peaks, despite good-faith efforts to 
push them together, it seems evident that the early- type Be stars cannot be all rotating with V eq = V CI i t . 



(15) 
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If the self-consistent grid of D values for the BO spectral type (shown in Figure 4a) was used instead of the 
lower-valued B9 grid, the modeled curves would be shifted further to the right in Figure 5 by about 18% 
(see vertical line in upper-right), thus making critical rotation even more of a mismatch. 

The quantitative degree of agreement between a simulated distribution $(vsini) (which is known ac- 
curately) and a coarsely sampled observed distribution is determined by using the x 2 diagnostic defined in 



where Nb = 11 bins, and where n b s and n mo d are the observed and modeled number counts in each bin, 
respectively. For the purposes of computing x 2 , the normalized distributions are multiplied by constant 
factors of either 92 or 93 in order for them to contain the same total number of stars as the observed 
subsample distributions. For the various models discussed below, Table 1 presents normalized ratios x 2 /V, 
where v = N\> — 1 is the effective number of degrees of freedom. Each resulting value of x 2 is related to the 
probability that the observed counts are drawn from the same "known" distribution from which the modeled 
counts are drawn. Assuming normally distributed uncertainties, this probability is given by 



where T(x) is the complete Gamma function. Naturally, when x 2 *C v the above probability approaches 
unity (i.e., the modeled distribution is a good match to the observations), and when x 2 » v the above 
probability is negligibly small. 

To locate the optimum range of V eq values for each subsample, the three free parameters of the trape- 
zoidal probability distribution ^(V^q) were varied until global minima in x 2 l v were found. This process was 
repeated 7 times, once for each of the following assumed values of the mean uncertainty level: — 0, 0.05, 
0.10, 0.15, 0.20, 0.25, and 0.30. In each case, the global minima in x 2 l v were located by: (1) automati- 
cally searching through the full three-dimensional parameter space, and (2) making contour plots of x 2 l v 
in several two-dimensional slices through the parameter space to make sure that no minima were missed. 
Table 1 lists the best-fit parameters and optimal values of x 2 /v 1 Q{x 2 \ v )> and for the two assumptions 
concerning V cr it discussed in § 2 (i.e., high-end and low-end limits). In all cases the best fits were found with 
either uq — 0.15 or 0.20. Lower values of tended to produce unrealistically sharp number distributions 
<b(vsmi), and higher values produced distributions that were too broad. This seems to be an independent 
verification that the observational uncertainties of the Yudin (2001) wsini values are about 10-20%. 

As seen in Table 1, there were two cases where no acceptable fits could be found to the observed number 
distributions: the L and VL subpopulations normalized by the low-end Chauville et al. (2001) critical rotation 
speeds. For these cases the peaks in the observed 3>(?;sini) distributions occurred at values of vsmi/V CI it 
larger than than 0.8, but the modeled statistical distributions — even for critical rotation — have peak values 
of v sini/V^rit no greater than ~0.7, as seen in Figure 5. If the Chauville et al. values of V C rit are correct, 
these late-type Be stars seem to be consistent with nearly critical rotation as well as possibly a weaker 
line-narrowing effect due to gravity darkening than has been modeled here. 

For ease of comparison with other analyses, Table 1 also gives the weighted mean ratios of equatorial 
rotation speed to critical rotation speed for each subpopulation, defined by 



§ 14.3 of Press et al. (1992): 




(19) 
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where x = V cq /V cr it is used for brevity. Taking into account that the five subpopulations have roughly 
equal numbers of stars, the mean ratio for the entire database of 462 stars is computed simply by averaging 
together the five values for each subpopulation. For the high-end (Cranmer) and low-end (Chauville et al.) 
choices of V cr i t values, the weighted mean values of (V cq /V cr it) are 0.684 and 0.854, respectively. The L 
and VL subpopulations that had no acceptable fits for ^(V eq ) in the low-end V CI it case were assumed to be 
completely critically rotating, so the above weighted mean value of 0.854 is an upper limit. Fremat et al. 
(2005) derived a most probable ratio of V cq /V a i t ~ 0.75 from an analysis of the 116 Be stars published by 
Chauville et al. (2001) plus 14 others. This value falls nearly halfway between the two mean ratios given 
above that were derived under the assumptions of lower and upper limiting cases for V^nt- 

Figure 6 presents a summary of the best-fit probability distributions 'J'(Veq) for each spectral- type sub- 
sample and for both the high-end and low-end calibrations of V a it- Figure 6 also compares the simulated and 
observed distributions of projected rotation speed $(wsini). Several features of these plots are noteworthy: 

1. There is a definite dependence of V m i n on spectral type. The hottest Be stars (subsamples VE and 
E) seem to have lower bounds on their rotation speed distributions that extend down to 40-50% of 
critical. This lower bound increases, as T e g decreases, to the point where the L and VL subsamples 
are consistent with nearly critical rotation. An extremely approximate fit to the dependence of V m i n 
on the stellar effective temperature is given by 

V m in/V crit « 0.28 tanh [13.3 - (T cff /1500 K)] + 0.72 (22) 

(see also Figure 8 below). 

2. The derived values of V max do not vary systematically with spectral type. Indeed, a decent broad-brush 
approximation would be to assume that Vmax ~ Krit, and that the rotation speeds of Be stars of a 
given spectral type range from a T e ff-dependent minimum value up to the critical rotation speed. 

3. The derived shapes of ^(V cq ) depend rather strongly on the adopted mean uncertainty level <7f. The 
higher the value of <t<; , the narrower the resulting best-fit probability distribution. This is understand- 
able because the simulated distributions <J>(usini) are modeled essentially as a convolution between 
two distributions: the intrinsic distribution of Veq values and the normally distributed spread of un- 
certainties (1 + (). If all of these distributions were Gaussian in shape, a specified width of $ could 
be obtained for an infinite number of choices for ^ and £, as long as the root-mean-squared sum of 
their widths equaled the desired width of <f>. Any future attempts to simulate these kinds of number 
distributions must be sure to model the observational uncertainties as accurately as possible. 

It is also worth noting that the choice of trapezoidal parameterization for ^(V eq ) (eq. [14]) was somewhat 
limiting because even the "best fitting" choices of the parameters did not yield perfect fits to the observed 
distributions $(wsini). For example, the strongly nonmonotonic behavior in the observed E subsample, 
for wsini/Vcrit ~ 0.6-0.7, limited the best-fit values of the probability Q{x 2 W) to be never larger than 
50%. A more complicated functional form for ^(y eq ) may have been able to fit this feature better, but the 
introduction of too many free parameters can lead to unphysically complex distributions. An attempt was 
made to vary the shape of \&(V C q) iteratively by using a "simulated annealing" algorithm — i.e., adopting 
randomized changes only if they resulted in lower values of \ 2 — following the spirit, if not the exact method, 
of Lucy (1974). The resulting best-fit distributions always tended toward sums of several (at least 3) sharply 
peaked subdistributions, with large ranges of intervening V eq having zero contribution. These distributions 
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were judged to be unphysical, and the more smoothly varying trapezoidal distribution was determined to be 
the best balance of realism and parameter flexibility. 

To develop a complete understanding of what ranges of V cq can actually be ruled out by this analysis, 
it is not enough to plot only the best-fit distributions. Figure 7 provides a summary of the goodness-of-fit 
probabilities Q(x 2 W) that were obtained both from the unconstrained parameter variation and also from 
other constrained searches of parameter space that examined only stars rotating faster than certain threshold 
values of V cq . Only results for the VE and M spectral subranges are shown, since the intermediate E subrange 
exhibits similar probability curves as the VE case and the L and VL cases are consistent with nearly critical 
rotation (and thus show no interesting behavior as the threshold lower limits are varied). These curves show 
how the probabilities for each subsample decrease when successively higher (i.e., more limited) ranges of 
Vmin and Vmax are allowed. The probabilities given in Table 1 for unconstrained searches of parameter space 
are shown in Figure 7 as the maximum values at the left edges of each plot. It is clear from Figure 7a that 
the VE number distribution <I>(vsini) cannot be fit adequately with values of Ve q that all exceed ~0.8 — no 
matter the choice of V CY it calibration. On the other hand, the necessity of subcritical rotation speeds for the 
M number distribution (in Figure 7b) depends sensitively on the adopted V cr it calibration. For the high-end 
(Cranmer) values, no good fits can be obtained when V cq is constrained to be greater than 0.8. For the 
low-end (Chauville et al.) values, though, critical rotation has just a slightly lower probability than nearly 
any degree of subcritical rotation and thus cannot be ruled out. 

The above constraints put on V m i n represent determinations of the threshold values of the equatorial 
rotation speed for the occurrence of the Be phenomenon. Figure 8 gives a summary of the T e g- dependence of 
the threshold rotation speeds from Table 1 and equation (22), as well as an indication of how V m i n changes 
for non-optimal values of (see caption for selection criteria). It should be emphasized that these values 
were all computed by assuming a relatively strong line-narrowing effect from von Zeipel gravity darkening. 
This assumption of D < 1 for fast rotation tends to transform a given distribution ^(V eq ) into a distribution 
$(i>sini) which is shifted more toward lower values of vsini than would have occurred if _£> = 1. If, though, 
D = 1 had been assumed in the above parameter optimization, the resulting values of Vmin would have all 
been smaller than the values given in Table 1 and Figures 6 and 8. Therefore, the values of Vmin derived here 
seem to be conservative upper bounds on the threshold for forming Be-star disks. If the gravity darkening 
effects were overestimated, then the Vmi n values may be lower, but it is difficult to imagine how they could 
be substantially higher. 

4. Departures from von Zeipel gravity darkening? 

The statistical results described in § 3 depend somewhat sensitively on the presumed form of gravity 
darkening used for the rapidly rotating stars. For stars with T e g > 8000 K it is traditional to assume that the 
mainly radiative energy transfer in the subsurface layers maintains the ideal von Zeipel (1924) linear scaling 
between emergent flux F and effective gravity g (i.e., F cx g). Cooler stars that exhibit subsurface convection 
are well-known to show weaker gravity darkening, with F cx g 4 P and the gravity darkening exponent (3 < 0.1 
(see, e.g., Lucy 1967; Osaki 1970; Anderson & Shu 1977; Claret 1998, 2000, 2003). Hot stars, though, may 
exhibit departures from the fully radiative gravity darkening value of /? = 0.25 if they are differentially 
rotating (Smith & Worley 1974), if they are in close binary systems undergoing mass transfer (Unno et 
al. 1994), or if there is some kind of shear-driven turbulence in their outermost layers that might help to 
redistribute the flux (e.g., Smith 1970; Smith & Roxburgh 1977). Some observational evidence exists from 
eclipsing binary light curves that /? could be greater than 0.25 in some cases (Kitamura & Nakamura 1988; 
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Nakamura & Kitamura 1992; Djurasevic et al. 2003). 

A general relationship between the bolometric flux F and the effective gravity g (the latter defined as 
the magnitude of the vector sum of gravity and the outward centrifugal force) is expressible as 



where 9 is the colatitude measured from the pole, <jb is the Stefan-Boltzmann constant, and [3 is the 
exponent defined in terms of effective temperature (i.e., T e ff oc g^). Small changes in as a function 
of the rotation rate of the star are neglected (see, however, Fremat et al. 2005). The parameter £4/3 is a 
normalizing constant that is a function of the rotation rate, but is not a function of 9. It is important to 
perform this normalization properly, since a key part of the v sin i analysis involves the comparison of stars 
modeled with and without gravity darkening. (In other words, it is important to know which latitudes are 
brighter, and which are dimmer, when compared with a star modeled without gravity darkening.) £4,3 is 
the surface-weighted integral of g 4 @ itself, and 



where 6 is the angle between the local effective gravity and the radius vector (cos<5 = —g r /g), and both g 
and the stellar radius i?* are functions of 9 (see, e.g., Slcttcbak 1949; Collins 1963; Maeder 1999). The two 
limiting cases of (3 = (no gravity darkening) and (3 = 0.25 (standard von Zeipel gravity darkening) give the 
surface area of the oblate star (S ) and the surface-weighted effective gravity (Si), respectively. For rigidly 
rotating Roche-model stars, Cranmer & Owocki (1995) presented a parameterized fit of Si as a function of 
the stellar rotation rate, and Cranmer (1996) gave a similar fit for So- For this paper, this function was 
also computed for an intermediate amount of gravity darkening between none and the standard amount 
{(3 — 0.125), and also for an extreme amount ({3 — 0.4) that was indicated by recent modeling of the light 
curves of early-type eclipsing binaries TT Aur and V Pup (Djurasevic et al. 2003). 

The line-narrowing factor D was computed for the two new choices for (3 and was compared with the 
values plotted in Figure 4 for the B9 spectral subtype. Figure 9a shows how the decrease in D with increasing 
V cq becomes stronger for larger values of (3. All curves are plotted for a constant inclination angle of i = 90°. 
Figure 9b illustrates how the simulated distribution of projected rotation speeds &j(vsini) shifts toward 
lower values for stronger gravity darkening exponents. This plot shows the distributions computed for a 
delta-function distribution of critical rotation speeds (Vj = V CT n) and thus is comparable to Figure 5. The 
curves in Figure 9b were computed assuming a mean uncertainty level of 0^ = 0.15. Note that for the case 
of no gravity darkening {(3 = 0) there would be a significant number of observations of v sin i values that 
exceed V a it- The fact that no such observations exist in the Yudin (2001) database implies that gravity 
darkening does need to be taken into account when processing the existing measurements of vsinz. Figure 
9b also shows the observed number distribution for the VL subsample of late-type Be stars, as normalized 
by the larger set of critical rotation speeds derived in § 2. The curve corresponding to (3 = 0.25 seems to 
be the most consistent with the observations, although values ±0.05 would also be reasonable. The shape 
of the observed probability distribution — especially at its upper end — may thus be a good diagnostic of the 
degree of gravity darkening that is present in a population of nearly critical rotators. 
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5. 



Linear Polarization 



Until this point there has not been much discussion of the dense circumstellar disks around Be stars 
that presumably exist only when V cq > V m [ n . A primary source of information about these disks is the 
measurement of linear polarization produced at broadband visible wavelengths (Coyne & McLean 1982; 
Bjorkman 2000b) by Thomson scattering of free electrons in the flattened envelope. The Yudin (2001) 
database contains 335 entries with nonzero values of both v sin i and the visible polarization fraction Py . 
Several aspects of the observed "triangular" distribution of Py versus v sin i were not easily explainable in 
terms of earlier models of Thomson scattering in circumstellar disks. Thus, the goal of this section is to use 
the stellar rotation properties derived above to produce a more accurate simulated distribution of polarization 
values and compare with the measured values. However, this is not an attempt to produce rigorous models 
of the physical properties of Be-star disks, but only a general consistency test for the validity of the derived 
ranges of rotation rate. 

Polarization values for a simulated distribution of Be stars were computed from a slightly modified form 
of the single-scattering approximation formulae of McDavid (2001). The required properties for each star 
were found by first computing three random quantities: (1) the inclination of the rotation axis, assumed to 
be identical to the disk axis, using the same procedure as outlined earlier; (2) the stellar spectral type, which 
was specified by a nonrotating value of T e g and was sampled from a cumulative probability distribution that 
was tabulated from the full Yudin (2001) database of 627 stars; and (3) the rotation rate V eq /V cr i t , which 
was sampled from a flat distribution (S — 0) between V m i n as given by equation (22) and V max — V cr i t . The 
remaining basic stellar parameters (R p , M*, and L„) were interpolated from the relationships derived in § 2, 
assuming that all stars are on the main sequence (luminosity class V). The equatorial values of the stellar 
radius i? cq and effective temperature T cq were computed assuming Roche oblatcness and ideal von Zeipel 
gravity darkening {[3 = 0.25). The disk temperature (Tdisk) was assumed to be constant and equal to 0.75 
times T cq (McDavid 2001). 

The disk itself was modeled as occupying a spherical sector around the equatorial plane with an opening 
half-angle a. The disk extends from its inner edge (assumed to be coincident with the star's equatorial 
radius R cq ) to infinity with a power-law dependence of electron density with radius, 



Following McDavid (2001), the constant value of rj = 3.1 is adopted for all stars. In an attempt to account 
for the large variation of fundamental parameters across the B-type spectral range, the inner disk density 
A^o was assumed to depend on T c g as 



This relation was derived from a linear fit (in logT c ff versus logTV^o) to a total of 23 measured inner disk 
densities for Be stars with spectral types ranging between B0. 5 and B8 (Waters et al. 1987; McDavid 2001). 2 
Figure 10 shows the T e g dependence of these values and compares the measurements to the above fit. The 
estimated mass densities p were computed under the assumption of complete hydrogen ionization; thus 
p = iVefflH where mn is the mass of a hydrogen atom. 



2 For another discussion of this trend, see Slcttcbak et al. (1992). Note, though, that van Kerkwijk et al. (1995) did not find 
such a trend in Be-star disk properties as derived from infrared observations. 



N e (r) = A e , (r/iZeq)"" . 
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To gauge the approximate validity of the derived inner disk densities, Figure 10 also plots two other 
densities that are expected to bracket these values from above and below. The photospheric mass densities 
are computed from the criterion that, in the stellar photosphere, the Rosseland mean optical depth should 
have a value of approximately one: 

t r w k rP H = 1 (27) 

where H is the photospheric scale height (proportional to T c $/g) and is the Rosseland mean opacity (in 
cm 2 g _1 ) interpolated from the extensive tabulation of Kurucz (1992). The resulting photospheric densities 
were compared with densities from detailed model atmospheres and were found to agree to within about 
±20% (R. L. Kurucz 2004, private communication). Figure 10 also shows an upper limit for the mass density 
at the sonic point of the equatorial outflow (i.e., the radius at which the radial flow speed equals the sound 
speed c s ), i.e., 

M 

Psonic — -. o ' 2 °) 

<\-KC s r l 

where the mass loss rate M was computed from the stellar wind parameterization of Vink et al. (2000). The 
use of these values presumes that the disk mass fluxes are of the same order of magnitude as the polar wind 
mass fluxes (see § 5.1 of Bjorkman 2000a). The radius r is set to R eq , though it should be noted that in a 
slowly expanding viscous decretion disk, the sonic radius may be much larger (e.g., Okazaki 2001) and thus 
the sonic density may be smaller. The sound speed c s is assumed to be constant inside the disk, 



c s = \Z-fkTdiskf ^rriH (29) 

where the adiabatic exponent 7 is 5/3, k is Boltzmann's constant, and fj, is the mean molecular weight of the 
gas (we assume \i — 0.6). The fact that the inner disk density falls between the photospheric density and 
the sonic density indicates that the base of the disk is most likely highly subsonic (i.e., close to the stellar 
surface), but is also sitting several scale heights above the photosphere. 

The opening half-angle a of the equatorial envelope was computed from the basic theory of Kcplerian 
accretion disks, with 

a = /tan-^c/Vfcep) (30) 
(see, e.g., Pringle 1981). In the above equation, Vk cp = (GM/r) 1 / 2 is the Keplerian azimuthal velocity at 
distance r, here computed at a fiducial distance of 2i? eq - The dimensionless constant / is an order-unity 
correction factor that is adjusted to produce an overall level of polarization similar to what is observed. The 
shape of the resulting distribution of Py versus v sinz does not depend strongly on the value chosen for /, as 
long as the same value is used for all stars in the simulated sample. In the models presented below, / = 3, 
which yielded a realistic range of a angles between 3° and 10°. The polarization from an axisymmetric 
disk is proportional to sin a cos 2 a, and for the small angles assumed here this is approximately a linearly 
increasing function of a. 

The above properties are the necessary inputs to the single-scattering formulae of McDavid (2001). The 
straightforward sin 2 i inclination dependence of that model, though, has been replaced by a more accurate 
relation that takes account of stellar occultation by a thin disk (i.e., eqs. [17] and [19] of Fox & Brown 
1991). The McDavid (2001) model was also simplified in two ways: (1) the relatively weak non-LTE effects 
for hydrogen were ignored, and (2) the polarization was computed only at the center of the V filter band 
(5500 A) instead of using the full bandpass function. The simulated values of v smi/V CI it for each star were 
computed by using equation (13) with <jq = 0.15. 

In Figure 11, the observed and simulated distributions of Py are compared with one another for sta- 
tistical samples of 335 stars each. The dotted triangular limits are shown only as a rough envelope of the 
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vsmi dependent spread of values, and are the same in both plots. Despite minor differences, the shapes 
of the observed and modeled distributions are similar. The rising upper envelope at low values of vsmi 
can be understood from the approximate sin 2 i dependence of the polarization, since the lowest projected 
rotation speeds correspond to the lowest inclination angles (see also McLean & Brown 1978). The decrease 
in the upper envelope for the largest values of wsini, however, is not so easy to interpret. Yudin (2001) 
suggested three possible explanations for this decrease in Py spread for the fastest rotators. Two of these 
effects (depolarization from the oblate gravity-darkened star itself, and lower intrinsic polarization for disks 
around the later-type stars that also have the highest values of V eq ) were not thought to be strong enough 
to produce the observed decreases. The third proposed effect was that the opening angle a may be inversely 
proportional to V eq . Whereas something like this may have been true for the wind-compressed disk model 
of Bjorkman & Cassinelli (1993), the thickness of a true viscous decretion disk seems most likely to become 
relatively uncoupled from the rotational properties of the central object (see eq. [30]). 

The simulated distribution of Py values shown in Figure lib does seem to match the observed strong 
decrease in spread at the largest values of vsmi in Figure 11a. This was found to occur because the models 
were computed using the gravity darkened equatorial temperature T oq as the argument of equation (26). 
The most rapidly rotating stars with the lowest values of T cq were thus given the lowest disk densities. A 
contributing factor was that these stars tended to be late-type Be stars with the lowest effective temperatures 
and the smallest radii. The stellar radius comes into play because the electron-scattering emissivity depends 
on the total number of electrons in the emitting volume (i.e., on the product N e R eq ). It remains to be 
shown whether the Be stars undergoing the most rapid rotation indeed have lower-density disks as predicted 
above. This seems somewhat counterintuitive, since whatever mechanisms that produce the disks are likely 
to grow in strength as the stellar rotation increases above the disk formation threshold (presumably at Vmin). 
However, for some combinations of stellar properties and rotation rate, the mass loss from the disk may come 
to dominate the angular momentum transport that feeds the disk, and lower densities may occur because of 
greater leakage to an equatorial wind. More work needs to be done to compare the derived disk properties 
of stars with similar spectral types but different rotation rates. 



6. Further Evidence for Subcritical Be-star Rotation 

The results obtained above imply that early-type Be stars can be rotating at significantly subcritical 
rates (down to roughly half of their critical rotation speed). This conclusion seems at odds with the recent 
statistical arguments of Townsend et al. (2004) that suggest nearly all Be stars could be rotating nearly 
critically. There are several additional pieces of evidence that can be brought to bear on this disagreement, 
and these ideas also imply subcritical rotation for at least some Be stars. 

Recent interferometric observations have put constraints on the inclination angles of some bright Be 
stars. Specifically, both Quirrenbach et al. (1997) and Tycner et al. (2004) measured an extreme axial ratio 
for the circumstellar envelope of C Tau that implies i > 70° . Using values of v sin i and V CI { t from Slettebak 
(1982) and Porter (1996), as well as assuming i = 70° in order to produce an upper limit, Tycner et al. 
(2005) estimated that V^ q /V^ r i t for this star is no more than ^0.52. Using the range of U cr it values computed 
in this paper for £ Tau and similarly dividing v smi/V CT it by sin 70° yields a slightly wider range of ratios: 
0.48 to 0.56. Assuming that gravity darkening was not taken into account in the derivation of wsini, the 
strongest possible correction would be to divide the above ratios by the smallest allowable value of D (i.e., 
0.68 for a B9 star), thus giving a range of V eq /V CI it for this star of 0.70-0.82. Given that the minimum value 
of D that corresponds to a Bl type star would be higher (^0.80), the subcritical nature of this star's rotation 
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appears definite. Future interferometric measurements are expected to bring about similar examples of firm 
subcritical rotation. 3 

In addition to rotationally broadened spectral lines, another potentially useful measure of near-critical 
rotation may be the large-scale shape of a star's spectral energy distribution (i.e., rotational "color effects"). 
It has been known for some time (e.g., Collins 1965; Kodaira & Hoekstra 1979) that gravity darkening can 
lead to more pronounced color variations in the ultraviolet than in the visible, though it is unclear how tightly 
these variations can be calibrated to provide a direct measurement of the rotation rate or the inclination. 
Stalio et al. (1987) examined far-ultraviolet flux distributions of several Be stars and was able to rule out 
extreme gravity darkening signatures. Stalio et al. made the preliminary conclusion that V eq /V CI it < 0.85 
for the observed stars. 

Collins & Sonneborn (1977) argued that critical rotation is unlikely for all Be stars because of the large 
predicted photometric shift above the main sequence. B-type stars rotating at 90% of their critical angular 
velocity (i.e., V cq /V cr it ~ 0.73) showed roughly a 1 magnitude increase in F-band absolute magnitude, while 
critical rotators showed more than a 2 magnitude increase. The presence of many Be stars at only ~1 
magnitude above the main sequence seemed to Collins & Sonneborn (1977) to be evidence against critical 
rotation. Note, though, that more recent calculations of photometric shifts due to rapid rotation have 
predicted a substantially weaker enhancement (see Figure 3 of Townsend et al. 2004). 

Mennickent et al. (1994) studied the detailed statistics of Ha spectral line shapes as a function of vsmi. 
From the observed spread of the measurements, they concluded that there is "...evidence for a considerable 
range of the true rotation velocities of Be stars: definitely there are intrinsically slow rotators among them." 

Finally, the observed variability of Be-star circumstellar envelopes may be related to the division found 
in this paper between subcritical rotation (for early-type Be stars) and nearly critical rotation (for late-type 
Be stars). Early- type Be stars are often observed to undergo strong outburst phases with an inferred rapid 
evacuation and refilling of the circumstellar disk region (see, e.g., Rivinius et al. 1998, 2001). Late-type Be 
stars are found to be comparatively quiescent, with the transitions between Be and normal-B phases being 
more gradual. If the rotation rates of the early-type Be stars are far below their respective V CT i t limits, 
they may require an impulsive mechanism of disk formation. Models of isotropic ejection from a point on 
the star — with some material being propelled forward into orbit and some propelled backward to fall back 
onto the star — may be needed to explain the disks around the hottest Be stars (Kroll & Hanuschik 1997; 
Owocki & Cranmer 2002). If, on the other hand, the coolest Be stars have nearly critical rotation rates, 
rapid outbursts as described above may not be needed for material to leak outwards into a viscous decretion 
disk (see also Clark et al. 2003; Owocki 2005). 

7. Summary and Discussion 

This paper contains a statistical analysis of the equatorial rotation rates of Be stars in Yudin's (2001) 
database. A new Monte Carlo forward modeling procedure was developed to simulate number distributions 



3 Several other stars discussed by Tycner et al. (2005) were found to have ratios Vi q /V cr it between 0.7 and 0.9, but applying 
the maximum possible correction for gravity darkening brings these stars to nearly critical rotation. The case of r\ Tau, for 
which Tycner et al. give a ratio of 0.53, seems to be anomalous because a main-sequence V cr it was used for this giant star; using 
the even the high-end value computed in this paper results in a ratio Vcq/V C rit of 0.80. This star is thus in the same category as 
the five others discussed by Tycner et al. (2005), excluding f Tau, and more detailed analysis must be performed to determine 
the star-specific corrections for gravity darkening. 
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of wsini for samples of strongly gravity-darkened hot stars. The parameters of the V cq distributions that 
best fit the observations were determined by a rigorous search of parameter space and a minimization of 
the appropriate \ 2 statistic. Although there were some variations in the resulting parameters depending on 
the assumed calibration of V CT n, the overall results seem robust. Early-type (07e-B2e) Be stars were found 
to exhibit a spread of equatorial rotation speed between a lower limit of 0.4-0.6 V cr it and an upper limit of 
critical rotation. Late-type (B3e-A0e) Be stars exhibit progressively narrower ranges of rotation speed as 
T e ff decreases; the lower limit rises gradually to 100% of critical rotation and the upper limit remains near 
the critical level. Uncertainties in V4 r it make it impossible to rule out critical rotation for Be stars later 
than about ~B3, but the result of substantially subcritical rotation for B0e-B2e stars appears to be firm. 
This paper (see Figure 8) thus seems to imply the existence of a rough dividing line in effective temperature 
between two qualitatively different types of Be phenomenon: 

1. For T e g > 21000 K the threshold ratio Vmin/V^rit above which a disk can form is well below unity. 
This seems to correspond to the regime of early-type Be stars that undergo violent (possibly pulsation- 
driven) outbursts that feed the circumstellar disk. 

2. For T c ff < 21000 if the threshold ratio Vmin/Krit ratio grows to unity as T e s decreases to the end of 
the B spectral range. The late-type Be stars at these temperatures are seen to be more quiescent than 
their early-type counterparts, and this may be due to the relative ease of generating a Keplerian disk 
for their closer proximity to critical rotation. 

The above results for early-type Be stars are in disagreement with other recent suggestions that the Be 
phenomenon is closely linked to critical or near-critical rotation (see, e.g., Maeder & Meynet 2000b). Pushing 
the V m - m threshold down to significantly subcritical rotation speeds makes it more difficult to explain the 
existence of Keplerian decretion disks around Be stars. 

The analysis presented in this paper can be improved in several ways to increase confidence in the 
results. The observational determination of wsini, V CT n, and possibly also the inclination angle i itself can 
be done more rigorously for an individual star when more than just one line profile is available (Stoeckley 
& Buscombe 1987; Reiners 2003; Fremat et al. 2005). For large wsini databases like that of Yudin (2001), 
the tabulated observational uncertainties should be utilized to a greater degree in producing comparable 
simulated samples. The assumption used in this paper of normally distributed uncertainties normalized by 
a single value of may have contributed to systematic biases that a more observationally guided procedure 
could correct. 

The various analyses performed here should also be done for catalogs of normal B-type stars that do 
not exhibit emission lines (see, e.g., Gkjbocki & Stawikowski 2000; Abt et al. 2002). It is still unclear if 
exceeding the minimum threshold rotation speed determined in this paper is both necessary and sufficient 
for the onset of the Be phenomenon; there may be many non-Be stars with V cq > V m i n . It is also not clear 
what sets the observed fraction of Be-star incidence as a function of the larger population of early-type stars 
(see, e.g., Zorec & Briot 1997; Penny et al. 2004). Firmer constraints on the rotation distributions of each 
population would be helpful. The definition of a truly "normal" B star is problematic, though, because a 
Be star can spend many years in a phase without a circumstellar disk and thus would have exhibited no 
Balmer emission over the entire time period of observations. There may be many early-type stars that will 
eventually become Be stars but are still classified as normal B stars. 

Our understanding of the rotation threshold for the onset of the Be phenomenon can be aided by incorpo- 
rating other kinds of observations that were not studied in this paper. Better interferometric measurements 
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that can constrain tightly both the stellar oblateness and the circumstellar disk geometry (McAlister et al. 
2005; Tycner et al. 2005) are becoming more widely available. In the long term, space-based missions like 
the proposed Stellar Imager may provide direct images of the oblate and gravity darkened stellar surfaces 
as well as firm constraints on differential rotation and macroturbulencc (Carpenter et al. 2004). Measuring 
the widths of lines that are formed in the narrow "boundary layer" above the subcritical photosphere, but 
below the inner edge of the Keplerian disk (e.g., Chen et al. 1989), may be a crucial probe of the physical 
origin of the Be phenomenon. 

I gratefully acknowledge Adriaan van Ballegooijen, Ian Howarth, Richard Townsend, and Stan Owocki 
for many valuable discussions, as well as Ruslan Yudin, Antonio Claret Dos Santos, and Bob Kurucz for 
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Administration (NASA) under grant NNG04GE77G to the Smithsonian Astrophysical Observatory. This 
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Table 1. Be Star Rotational Properties by Subpopulation 



VE 



E 



M 



VL 



Definitions: 

min (T cff ) (K) 24000 20500 

max (T cff ) (K) 36000 24000 

approx. spectral types 07-B1.5 B1.5-B2.5 

number of stars 93 93 

Observed statistics (Cranmer's Vcrit): 

mean (v sin i/V crit ) 0.4069 0.4335 

std. dev. (vami/V crit ) 0.1527 0.1401 

skewness (v sin i/Vcrit) 0.2625 0.3260 

kurtosis (vsmi/V cli t) 0.1242 -0.4883 

Observed statistics (Chauville et al. V^rit): 



18500 
20500 
B2.5-B3.5 
92 



0.4740 
0.1693 
-0.1048 
-0.1489 



13200 
18500 
B3.5-B6 
92 



0.5269 
0.1619 
-0.4354 
-0.5136 



10000 
13200 
B6-A0 
92 



0.5863 
0.1490 
-0.1232 
-0.3359 



mean (v sinz/Vcrit) 


0.4926 


0.5225 


0.5812 


0.6795 


0.8957 


std. dev. (t> sin i/Vcrit) 


0.1837 


0.1661 


0.2061 


0.2144 


0.2544 


skewness (v sin i/V cr it) 


0.2230 


0.1964 


-0.4215 


-0.7592 


-1.2430 


kurtosis (vsini/Vcrit) 


0.3961 


-0.7920 


-0.5134 


-0.3913 


-0.4913 


Best-fit V eq distributions 


(Cranmer's Vc r i t ): 








Vmin/Vcrit 


0.4150 


0.3423 


0.5050 


0.7219 


0.9084 


^max/ Vcrit 


0.7862 


0.7955 


0.9118 


0.7454 


0.9121 


s 


-0.9500 


0.0886 


-0.7645 


0.6904 





min (x 2 /V) 


0.39 


1.31 


0.58 


0.44 


0.57 


Q(x»,in% 


95.3 


21.6 


83.5 


92.5 
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Fig. 1. — Critical rotation speed as a function of spectral type for stars of luminosity class III (dotted line), 
IV (dashed line), and V (solid line). Symbols denote values of V a n taken from other papers (see labels 
above). 



-27- 



3 
2 
1 



(a) all 



0.0 



0.2 



0.4 



0.6 



0.8 



1.0 



2 
1 



(b) VE 



0.0 



c 

T3 



(c) E 



>, 0.0 

OS 
-Q 
O 

& 



(d) M 



0.0 



0.2 



0.2 



0.2 



0.4 



0.4 



0.4 



0.6 



0.6 



0.6 



0.8 



0.8 



0.8 



1.0 



1.0 



1.0 



2 
1 



(e) L 



0.0 



0.2 



0.4 



0.6 



0.8 



1.0 



2 
1 



(f) VL 



0.0 0.2 0.4 0.6 

(v sin i)/V crit 



0.8 



1.0 



Fig. 2. — Observed number distributions of vsmi/V a i t for: (a) the full Yudin (2001) database of 462 stars 
having nonzero v sini; (b)-(f) the five subpopulations, defined by the T e g bounds given in in Table 1. These 
ratios are computed using the high-end (Cranmer) calibration for V^ r it . Also shown are the means (dashed 
lines) and ±1<t standard deviations (dotted lines) for each of the independently normalized distributions. 
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Fig. 3.— Schematic continuum intensity images of a B9 V star rotating rigidly with V cq = 0.95 V cr it and 
viewed from an inclination angle i = 75° measured from the pole. Both images are derived with the same 
linear limb darkening, but they differ in that (a) uses the centrifugal force to alter the surface equipotential 
surfaces (Roche oblateness) and redistributes the radiative flux in proportion to the effective gravity (von 
Zeipel gravity darkening), and (b) does not. Panel (c) shows resulting absorption lines computed as described 
in the text: the solid curve corresponds to image (a); the dashed curve corresponds to (b). The central depth 
of profile (b) was renormalized to match that of (a) to more clearly compare the line shapes. 
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Fig. 4. — Contour plots of the line narrowing ratio D due to gravity darkening, for main sequence BO (a) 
and B9 (b) models. The values of the constant levels are labeled above (white denotes D~l). 
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Fig. 5. — Statistical distributions <!>.,■ (u sin i) computed under the assumption that all stars are rotating 
at V^ r it and using the strongest (i.e., B9) gravity darkening grid for D. The simulated observations were 
made with varying levels of random uncertainty: er^ = (solid line), 0.05 (dashed line), 0.1 (dot-dashed 
line), 0.2 (dash-triple- dot line), and 0.3 (dotted line). Also plotted (gray bars) is the observed distribution of 
vsmi/V CY it values for the 93 stars in the VE subpopulation, normalized using the fit to the Chauville et al. 
(2001) low-end V CT it values. A vertical line in the upper-right shows how the distributions would be shifted 
to higher values of vsini if the consistent gravity darkening grid for the VE stars (i.e., B0) would have been 
used. 
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Fig. 6. — Comparisons of observed (solid lines) and simulated (dashed lines) probability distributions of the 
projected rotation speed &(v sini), as well as the best-fit distributions of intrinsic rotation speed ^(V eq ) 
(dotted lines) for the 5 spectral subpopulations defined in Table 1. (a)-(e): calculations assuming the high- 
end (Cranmer) calibration for V cr i t derived in § 2. (f)-(h): calculations assuming the low-end (Chauville et 
al.) calibration for V a a- The L and VL models for the latter calibration did not yield acceptable solutions 
and are not shown. 
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Fig. 7. — Threshold goodness-of-fit probabilities Q for the (a) VE and (b) M spectral subranges, computed 
for a range of constrained searches of parameter space. Each data point represents the maximum probability 
allowed when V m i n and V max were required to exceed the threshold values plotted as the abscissa. Parameter 
searches were performed assuming both the high-end (Cranmer) V a n calibration (solid lines) and the low-end 
(Chauville et al.) V cr i t calibration (dashed lines), both with kept constant at 0.20. 
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Fig. 8. — Minimum allowed equatorial rotation speeds for the observed populations of Be stars, plotted as 
dimcnsionless ratios V m i n /V cr it versus T e g . The empirical limits determined from x 2 minimization (see Table 
1) are plotted as horizontal bars corresponding to the high-end Cranmer V cr n calibration (solid lines) and 
the low-end Chauville et al. V^rit calibration (dotted lines). The gray bars show the range of V min values 
that occur when is varied between and 0.30 in steps of 0.05, and only high-probability values with 
Q(x 2 W) > 10% are kept. Also shown is the approximate fitting formula given by equation 22 (dashed line). 
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Fig. 9. — (a) Line narrowing ratio D versus equatorial rotation speed for a constant inclination angle of 
90° and a range of gravity darkening exponents f3 (see labels above) . (b) Statistical distributions $j (v sin i) 
computed for the same range of (3 values as in panel (a), assuming that all stars arc rotating at V CT it and 
have a mean uncertainty level <jq = 0.15. Also plotted {gray bars) is the observed distribution of v smi/V CT it 
values for the 92 stars in the VL subpopulation normalized by the high-end Cranmer calibration for V^nt 
(see also Fig. 2f). 
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Fig. 10. — Mass densities (in g cm -3 ) computed for main sequence nonrotating B stars at various heights, 
versus stellar effective temperature: optical-depth-unity photosphere (solid line), an upper limit estimate for 
the equatorial sonic point of disk outflow (dashed line), and the fitting formula for the inner disk density 
given in equation (26) (dotted line). Also shown are 8 measurements of inner disk densities from Table 2 of 
McDavid 2001 (triangles) and 15 measurements from Tables 3 and 5a of Waters et al. 1987 (squares). 
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Fig. 11. — Correlation between visible-band polarization Py and projected rotation speed vsini, shown for: 
(a) 335 observed stars with nonzero measurements for both values from the Yudin (2001) database, and (b) 
335 simulated stars, using a modified version of the method of McDavid (2001). Rough upper bounds to the 
spread of Py values are also plotted (dotted lines). 



